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Abstract 

In reality, many real-world networks interact with and depend on other networks. We develop an analytical framework for 
studying interacting networks and present an exact percolation law for a network of n interdependent networks (NON). We 
present a general framework to study the dynamics of the cascading failures process at each step caused by an initial failure 
occurring in the NON system. We study and compare both n coupled Erdos-Renyi (ER) graphs and n coupled random regular 
(RR) graphs. We found recently [Gao et. al. arXive: 1010.5829] that for an NON composed of n ER networks each of average 
degree k, the giant component, P x , is given by = p[l — exp(— kP^)]" where 1 — p is the initial fraction of removed nodes. 
Our general result coincides for n = 1 with the known Erdos-Renyi second-order phase transition at a threshold, p — p c , 
for a single network. For n = 2 the general result for corresponds to the n = 2 result [Buldyrev et. al., Nature , 464, 
(2010)]. Here we show for an NON composed of n coupled RR networks each of degree k, that the giant components is given by 

Poo = p{l - {p 1/n P£~ 1)/n {(l - (P 00 /p) 1/ ") (fe ~ 1)/fe - 1] + l} fe }™. Similar to the ER NON, for n = 1 the percolation transition 
at p c , is of second order while for any n > 1 it is of first order. The first order percolation transition in both ER and RR (for 
n > 1) is accompanied by cascading failures between the networks due to their interdependencies. However, we find that the 
robustness of n coupled RR networks of degree k is dramatically higher compared to the n coupled ER networks of average 
degree k. While for ER NON there exists a critical minimum average degree k — k m i u , that increases with n, below which the 
system collapses, there is no such analogous fc m i n for RR NON system. For any k > 2, the RR NON is stable, i.e., p c < 1. This 
is due to the critical role played by singly connected nodes which exist in ER NON and enhance the cascading failures but do 
not exist in the RR NON system. 



I. INTRODUCTION 

In our modern world, infrastructures, which affect all 
areas of daily life, are usually interdependent. Exam- 
ples include, electric power, natural gas and petroleum 
production and distribution, telecommunications, trans- 
portation, water supply, banking and finance, emergency 
and government services, agriculture, and other funda- 
mental systems and services that are critical to the se- 
curity, economic prosperity, and social systems, shown 
in Fig.l. Although urban societies rely on each of the 
individual infrastructures, recent disasters ranging from 
hurricanes to large-scale power blackout and terrorist at- 
tacks have shown that significant dangerous vulnerabil- 
ity is due to the many interdependencies across differ- 
ent infrastructures [1-5]. Infrastructures are frequently 
connected at multiple points through a wide variety of 
mechanisms, such that a bidirectional relationship exists 
between the states of any given pair of networks, as shown 
in Fig. 1 and 2 [1]. For example, in California, electric 
power disruptions in early 2001 affected oil and natural 
gas production, refinery operations, pipeline transport of 
gasoline and jet fuel within California and its neighbor- 
ing states, and the movement of water from northern to 
central and southern regions of the state for crop irriga- 
tion. Another dramatic real-world example of a cascade 
of failures is the electrical blackout that affected much 
of Italy on 28 September 2003: the shutdown of power 
stations directly led to the failure of nodes in the Su- 
pervisory Control and Data Acquisition (SCADA) com- 



munication network, which in turn caused further break- 
down of power stations [5, 6]. Identifying, understand- 
ing, and analyzing such interdependencies are significant 
challenges. These challenges are greatly magnified by the 
breadth and complexity of our modern critical national 
interdependent infrastructures [4]. 

In recent years we observed important advances in the 
field of complex networks [7-19]. The internet, airline 
routes, and electric power grids are all examples of net- 
works whose function relies crucially on the connectivity 
between the network components. An important prop- 
erty of such systems is their robustness to node failures. 
Almost all research has been concentrated on the case 
of a single or isolated network which does not interact 
with or depend on other networks. Recently, based on 
the motivation that modern infrastructures are becom- 
ing significantly more dependent on each other, a system 
of two coupled interdependent networks has been stud- 
ied [6, 20, 21]. A fundamental property of interdependent 
networks is that when nodes in one network fail, they may 
lead to the failure of dependent nodes in other networks 
which may cause further damage in the first network and 
so on, leading to a global cascade of failures. Buldyrev et 
al. [6] developed a framework for analyzing the robustness 
of two interacting networks subject to such cascading fail- 
ures. They found that interdependent networks behave 
very different from single networks and become signifi- 
cantly more vulnerable compared to their noninteracting 
counterparts. 

In many realistic examples, more than two networks 
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depend on each other. For example, diverse infrastruc- 
tures, such as water and food supply, communications, 
fuel, financial transactions, and power stations are cou- 
pled together [1, 3, 5, 22]. Understanding the vulnerabil- 
ity due to such interdependencies is a major challenges 
for designing resilient infrastructures. 



FIG. 1: Illustration of the interdependent relationship 
among different infrastructures [1]. These complex rela- 
tionships are characterized by multiple connections be- 
tween infrastructures, feedback and feedforward paths, 
and intricate, branching topologies. The connections 
create an intricate web that, depending on the charac- 
teristics of its linkages, can transmit shocks throughout 
broad swaths of an economy and across multiple infras- 
tructures. It is clearly impossible to adequately analyze 
or understand the behavior of a given infrastructure in 
isolation from the environment or other infrastructures. 
Rather, one must consider multiple interconnected in- 
frastructures and their interdependencies. For example, 
the reliable operation of modern infrastructures depends 
on computerized control systems, from SCAD A systems 
that control electric power grids to computerized systems 
that manage the flow of railcars and goods in the rail in- 
dustry. In these cases, the infrastructures require infor- 
mation transmitted and delivered by the communication 
infrastructure [1]. 



of n coupled randomly connected networks. Our result 
generalizes the known Erdos-Renyi (ER) [24-26] result 
as well as the random regular (RR) result for the giant 
component of a single network, and shows that while for 
n = 1 the percolation transition is a second order, for 
n > 1 cascading failures occur and the transition be- 
comes a first order transition. Our results for n inter- 
dependent networks show that the classical percolation 
theory extensively studied in physics and mathematics is 
in fact a limited case of the rich, general, and very dif- 
ferent percolation law which exists in realistic interacting 
networks. 

Additionally, we find: 

(i) for any loopless topology of NON, the critical per- 
colation threshold and the giant component depend only 
on the number of networks involved and their degree dis- 
tributions but not on the inter-linked topology (Fig. 2), 

(ii) the robustness of NON significantly decreases with 
n, and 

(iii) for a network of n ER networks all with the same 
average degree k, there exists a minimum degree £; m j n (n) 
increasing with n, below which p c = 1, i.e., for k < k mm 
the NON will collapse once any finite number of nodes 
fail. The analytical expression for k m i n (n) generalizes the 



known result k„ 
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1 for ER below which the network 



collapses. In sharp contrast a NON composed of RR 
networks is significantly more robust. In the RR NON 
case there is no k m - m below which the NON collapses. 
This is due to the multiple links of each node in the RR 
system compared to the existence of singely connected 
nodes in the ER case. We also discuss the critical effect 
of singly connected nodes on the vulnerability of the NON 
ER structure. 




line-like NON star-like NON tree-like NON 



FIG. 2: (color online) Three types of loopless NONs com- 
posed of five coupled networks all have same percolation 
threshold and same giant percolation component. 



We study here a model system [23], comprising a net- 
work of n coupled networks, where each network con- 
sists of N nodes (See Fig. 2). The N nodes in each 
network are connected to nodes in neighboring networks 
by bidirectional dependency links, thereby establishing 
a one-to-one correspondence. We apply a mathematical 
framework [23] to study the robustness of tree-like "net- 
work of networks" (NON) by studying the dynamically 
process of the cascading failures. We find an exact an- 
alytical law for percolation of a NON system composed 



II. THE DYNAMIC PROCESS OF CASCADING 
FAILURES 

To model an interdependent NON we consider for sim- 
plicity and without loss of generality, n networks each 
having N nodes. We study the percolation of n networks 
connected in a loopless structure, the structure of the 
NON can be, e.g., a line, a star or a tree as shown in 
Fig. 2. Each node in Fig. 2 represents a network, and 
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each link between two networks i and j denotes the ex- 
istence of a one-to-one dependencies between the nodes 
of the linked networks. The functioning of one node in 
network i depends on the functioning of one and only 
one node in network j G {1, 2, n}, i ^ j), and 
vice versa (bidirectional links). We assume that within 
network i, the nodes are randomly connected by A^-links 
with degree distribution P(fci), where ki is the average 
degree of network i. 

The root of the NON is the network from which frac- 
tion 1 — p of nodes are removed due to random failure. 
Before showing the dynamic of the cascading failures, we 
present the following three definitions, (i) We define the 
distance matrix Dij as the distance form network i to 
network j in the NON. (ii) Shell j is a set Lj whose net- 
works are at distance j from the root network, where 
j G [0, s] and s is the total number of shells. In the fol- 
lowing example, we use ij to denote network i in shell j, 
e.g., ij G Lj. Note that in shell 0, there is one and only 
one network 1 . (Hi) gi(x) is the generating function of 
network i [27], which reflects the topology of network i 
and satisfies 

gi{x) = 1 - G 0;i (x, fi). (1) 
where Go,i{x,fi) satisfies [6, 27, 28] 

oo 

Go,i(z, fi) = G ,i(xfi + l-x) = J2 Pi(k)(xfi + 1 - x) k , 

(2) 

and 

G 1 (x,f) = G' (x,f)/G' (l) = f. (3) 

Next, we show analytically the steps in the dynamics 
of the cascading failures as demonstrated in Fig. 3. 




(d) t=3,5,7,9,... (e) t=4,6,8,10,... 

FIG. 3: The dynamic of cascading failures. In this figure, 
each node represents a network. The arrow (on the link) 
illustrates the damage spreading from the root network 
to the whole NON shell by shell. 



Step (0): At t = (Fig. 3(a)), we begin by randomly 
removing a fraction 1 — p of nodes from the root network 
(network lo), and removing all the Ai -links connected 
to these removed nodes. Next we remove all the nodes 
that become disconnected to the largest component of 
network lo- Thus at t = 0, for the network lo, the frac- 
tion of remaining nodes in network lo after the initial 
failure is xq.i = P an d the fraction nodes in the giant 
component // ,i = pgi (p). 

Step (1): At t = 1, the root network spreads its dam- 
ages to all its neighboring networks i\ G L\ (Fig. 3(b)). 
So we remove all nodes in networks i\ that are connected 
to the removed nodes in network lo and then remove all 
the nodes not in the giant components of networks i\. At 
t = 1, the failure of networks i\ is equivalent to a random 
removal of the fraction of \—x\^ x nodes from networks i\ 
[6], where Xi ^ = pgi (xo i i ), and the giant component 
of network i x is = x lt i 1 g il (xi t i 1 ). 

Step (2): At t = 2, the networks i\ € L\ reflects their 
damages back to the root network and spreads their dam- 
ages to all their neighboring networks j2 G ^2 (See Fig. 
3(c)). So we remove all nodes in networks lo and ji that 
are connected to the removed nodes in networks i\ and 
then removing all the nodes not in the giant components 
of networks i±. Again the failure of network lo is equiva- 
lent to a random removal of the fraction of 1 — X2,i nodes 
from networks 1 , where x 2 ,i = pl\ iieLl 5u(a;i,u); the 
failure of networks j2 is equivalent to a random removal 
of the fraction of 1 — x 2 j 2 nodes from networks j 2 , where 
x z,j2 = Pgio( x o,io)9ii( x i,ii) f° r networks i\ that are 
linked to networks j 2 ■ 

Step (3): At t = 3, the root network spreads its fur- 
ther damages to the networks i\ in shell 1 again, the 
networks j 2 in shell 2 reflect their damages back to the 
neighboring networks i\ in shell 1, and to the neigh- 
boring networks 113 in shell 3 as shown in Fig. 3(d). 
A network in i\ receives the damages information x 2i \ 
from network lo, x 2 .i 2 from networks j 2 that are linked 
to networks i\, and x\. Vl from networks v\ in shell 1 
where the networks v\ are the neighboring networks of 
i\S neighboring networks, i.e., the distance between net- 
works i\ and networks v\ is 2. Thus we can obtain that 
x 3 ,h = P9i {x2 } i )U.j 2 9h{x2,j 2 )Uv 1 9v 1 (xi, Vl ). Simi- 
larly, we obtain the failure of networks U3 G £3 to be 
^3,« 3 = Pgi (x 0tla )g ll (x 1 ^ 1 )g j2 (x 2 j 2 ) 1 where networks u 3 
in shell 3 are connected to networks j 2 in shell 2 and 
networks j 2 are connected to networks i\ in shell 1. 

We continue the cascading process step by step [see 
Figs. 3(d) and (e)] until the convergence step, t = r, 
when no further nodes and links removal occurs. Accord- 
ingly, we investigate the dynamically cascading process 
of our model of the loopless NON. First we initialize the 
NON as 

^o,io =P)Mo,i =P9i (p), (4) 

and 

gyixt,^) = MM, =l,t<j. (5) 
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FIG. 4: How does the damage spread in a NON system? 
In this figure, each node represents a network. When 
looking at network 12 for example, it becomes damaged 
at t = 2k + l (k = 1, 2, 3, ...). It receives the damage from 
network 8 at t = 2k+3, because network 8 gets damage at 
t = 2 for the first time and its damage spreads to network 
12 when t = 5 for the first time, which agrees with Eqs. 
(4-8) that network i receives damage from network j if 
and only if t — Dij > D\j . 



Furthermore, when t — > oo, ^ = /i^ and xt-D i:j ,i = 
x t .i = Xi. So from Eqs. (4-8) we can obtain that 



x%=V 9j(xj), (9) 

and 

n 

Moo =p*[[gj(xj). (io) 

3=1 

We also demonstrate how does the damage spread in 
another example of NON shown in Fig. 4. 

In Figs. 5 and 6 we compare our theoretical results 
Eqs. (9) and (10) with simulation results for 3 different 
types of networks, ER networks, RR networks and SF 
networks. We find that while the dynamics is different for 
the three topologies shown in Fig. 2, the final = jtioo.i 
is the same as predicted by the theoretical results, Eq. 
(27). 

Next, we study the case for n coupled networks where 
all networks are with the same degree distribution spec- 
ified by the generating functions Gq(x, /) (Eq. (2)) and 
Gi(x, f) (Eq. 3)). By substituting Eq. (1) and (2) into 
the Eq. (9) and introducing the parameter z = xf+l—x, 
we obtain 



x - p(l - Go(z))"- 1 (11) 

From Eqs. (3), (11) and (10), the equations for mutual 
giant component become 



Thus we can obtain that the giant component of net- 
work ij in shell j, fj^^ at step t satisfies 

= Xt^gijix^J € C(t), (6) 



Mm, = Mt-M^j i C{t). (7) 

where C(t) satisfies the sequence 

r c(o) = {o} 

C(1) = {1} 
C(2) = {0,2} 
C(3) = {1,3} 
C(4) = {0,2,4} 
C(5) = {1,3,5} 
C(6) = {0,2,4,6} 
(7(7) = {1,3, 5, 7} 
(7(8) = {0,2,4,6,8} 



and fl^i satisfies 



n 

x t ,i=p JJ gj{x t - Dij ,i}- (8) 

3=1 



i = (l-Go^r-Hi-Gi^)) ^ (12) 

p 1—2! 

and 

fio=P(l-G°(^= (1 ~f^ ) g ) " g) - (13) 



III. THE CASE OF NON COMPOSED OF n ER 
NETWORKS 

The case of NON of n Erdos-Renyi (ER) [24-26] 
networks with average degrees k\, k 2 , ...fcj, k n can be 
solved explicitly [23]. In this case, the generating func- 
tions of the n networks are [28]. 

G h i(x) = G Q ,i(x) = exp[fci(a; - 1)]. (14) 

Accordingly, we obtain that the generating function 
gi(xi) satisfies 

gi{xi) = 1 - exp[fei^(/i - 1)], (15) 
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theory of star like NON 
theory of line like NON 
theory of tree like NON 
simulation of star like NON 
simulation of line like NON 
simulation of tree like NON 

(a) 




FIG. 5: (a) Simulation results of the giant component 
of the root network ^ t ,i after t cascading failures for 
three types of NON composed of 5 ER networks shown 
in Fig. 2. For each network in the NON, N = 100, 000 
and k = 5. The value of p chosen is p — 0.85, and the 
predicated threshold p c = 0.76449 (from Eqs. (22) and 
(24)). All points are the results of averaging over 40 re- 
alizations. Note that while the dynamics is different for 
the three topologies, the final Poo = /ioo.i is the same, 
i.e., the final does not depend on the topology of the 
NON. (b) Simulations of the giant component, [M t i, for 
the tree-like NON [Fig. 2] with the same parameters as 
in (a) but for p = 0.755 < p c = 0.76449. The figure 
shows 50 simulated realizations of the giant component 
left after t stages of the cascading failures compared with 
the theoretical prediction of Eqs. (4)-(8). 
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FIG. 6: (a) Simulation results for the giant component of 
the root network fj, tt \ after t cascading failures for three 
types of NON composed of 5 RR networks. The struc- 
tures of the NON are as shown in Fig. 2. For each net- 
work in the NON, N = 100, 000 and k = 5 . The value 
of p chosen is p — 0.65, and the predicated threshold 
p c = 0.6047 [from Eqs. (41) and (43)]. The points are 
the results of averaging over 40 realizations. It is seen 
that while the dynamics is different for the three topolo- 
gies, the final Poo = /ioo,i is the same, i.e., the final Poo 
does not depend on the topology of the NON. (b) Simu- 
lation results of the giant component of the root network 
Ht,i after t cascading failures for three types of NON 
composed of 5 SF networks shown in Fig. 2. For each 
network in the NON, N = 100, 000 and A = 2.3 . The 
value of p chosen is p = 0.85(above p c ). The points are 
the results of averaging over 40 realizations. We also can 
see here that while the dynamics is different for the three 
topologies, the final P>o = Hoo,i is the same. 



where fi = exp[kiXi(fi - 1)] and thus gi(xi) = 1 - 
Using Eq. (9) for Xi we get 

n 

fi = eM-Ph l[a - /,•)], i = l,2,...,n. (16) 

By introducing a new variable r = f^ kt , i = 1,2, n 
into Eq. (16), we can reduce the n equations to a single 
equation, 



r = cxp[— pJ^J(l — r k 



(17) 



which can be solved graphically for any p. For small p, 
Eq. (17) has only the trivial solution r = 1. This case 
corresponds to the absence of the mutual giant compo- 
nent and hence to the complete fragmentation of the net- 
works. As p increases a nontrivial solution r < 1 emerges 
at some critical value of p = p c . The critical case corre- 
sponds to the tangential condition: 



1 = J- exphpjja- »■*')]■ 



(18) 



i=i 



Thus, the critical value of r satisfies a transcendental 
equation 



ft 

n(l-r*) 

[ i=i 
r = cxp | - — 



i=l 



n (i-r*or 



(19) 



From Eqs. (16) and (18) we can obtain the critical 
percolation threshold p c and the the value of /Zoo at p c as 



n n _ ^ 

1=1 j'=i,j¥» 



and 



(20) 



(21) 



i=l 



If p < Pc, Eqs. (16) have only the trivial solutions 
(fi — 1) and /loo = 0. When the n networks have the 
same average degree k, ki — k (i = 1, 2, n), we obtain 
from Eq. (16) that f c = fi(p c ) satisfies [23] 



fc 



(22) 



This solution can be expressed in terms of the Lambert 
function W(x) [31, 32], 



f c = -[nW(~e-^]-\ 



(23) 



Once f c is known, we obtain p c and /ioo.„ = Poo at p c by 
substituting ki — k into Eqs. (20) and (21) 
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p c = [nfc/ c (l-/ c )("- 1 )]- 1 



(24) 



fcminM = [n/c(l - /c) (n - 1} ] 



(26) 



and 



P — 



l-/c 

nkf c 



(25) 



For n = 1 we obtain the known results p c = 1/fc and 
P<x = at p c (representing the second order transition) 
of Erdos-Renyi [24-26]. Substituting n = 2 in Eqs. (24) 
and (25) one obtains the exact results derived in [6]. Note 
that for all n > 1 we obtain > at p c representing a 
first order nature of the percolation transition. For the 
behavior of p c [Eq.(24)] for large n see Appendix A. 



Note that Eq. (26) together with Eq. (22) yield the 
value of £; m i n (l) = 1 for n = 1, reproducing the known 
ER result, that (fc) = 1 is the minimum average degree 
needed to have a giant component. For n = 2, Eq. (26) 
yields the result obtained in [6], i.e., fc m ; n = 2.4554. 

When the n networks have the same average degree 
k, ki — k (i — 1, 2, n), using Eqs. (17) and (21) we 
obtain the percolation law for the order parameter, the 
size of the mutual giant components for all p values and 
for all k and n [23], 



/Jo 



p[l - exp(-kP a 



(27) 



A. The minimum degree k and the giant compo- 
nent Poo(p) 



The solutions of equation (27) for several n-values are 
shown in Fig. 8(a). Results are in excellent agreement 
with simulations. The special case n = 1 is the known ER 
second order percolation law for a single network [24-26] . 




n 



FIG. 7: The critical fraction p c for different k and n 
for (a) ER NON system and (b) RR NON system. The 
results for ER NON system are obtained by Eqs. (22) and 
(24), while the results of RR NON system are obtained 
from the solution of Eqs. (41) and (43). The results are 
in good agreement with simulations. In the simulations 
p c was calculated from the number of cascading failures 
which diverge at p c [33] 

Here we show that while for the ER NON system there 
exists a minimum degree fc m ; n (n) below which the NON 
collapses, in the RR NON such a fc m ; n does not exists and 
the RR NON is stable (p c (n) < 1) for any n. Thus, the 
RR NON is significantly more robust compared to the ER 
NON due to the critical role of the singly connected nodes 
on the vulnerability of the NON system. To analyze p c as 
a function of n, we find f c from Eq. (22) and substitute 
it into Eq. (24), and we obtain p c as a function of n for 
different k values, as shown in Fig. 7(a) for the ER NON 
system. It is seen that the NON becomes more vulnerable 
with increasing n or decreasing k (p c increases when n 
increases or k decreases). Furthermore, for a fixed n, 
when k is smaller than a critical number fc m »„(n), p c > 1 
meaning that for k < fe m » n (n), the NON will collapse 
even if a single node fails. From Eq. (24) by substituting 
p c = 1, we get the minimum of fc as a function of n 



B. The case of NON with different average de- 
grees 

Next, we study the case where the average degrees of 
all n networks is not the same. Without loss of generality 
we assume that m networks have the same average degree 
(fc)2, and other n — m networks have the same average 
degree (k)\. We define a — (fc)i/(fc)2 where < a < 1. 
Using Eqs. (16-19) we can show that f c satisfies 

/ = ex P [ (/c-D(l-/c 1/a ) , 

(n - m)/ c (l - fc /a ) + mfl ,a {\ - f c )/a 

(28) 

Results for p c and the mutual giant component for dif- 
ferent values of m, and alpha are shown in Fig. 11. The 
case of (fc)i <C (A;) 2 is interesting, since in this limit the 
m networks with large (fc)2 due to their good connectiv- 
ity can not cause further damage to the n — m networks 
with (fc)i. Thus the NONs system can be regarded as 
only n — m networks. Indeed, when a — >• 0, f c satisfies 

f c = eT^k. (29) 
And then equation of p c and obtained as 

p c = [(n - m)(fc) 1 / c (l - /e)^-™" 1 )]- 1 , (30) 

and 

Moo = 7 T77T 7~- (31) 

(n - m)(fc)i/ c 

Equations (29-31) are indeed the same as Eqs. (22, 24, 
and 25) where n is replaced byn-m. This result is seen 
also in Fig. 11, where the limit of a — yield the same 
results as for a = 1 for n — m networks. 
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When (fc)i <C {k) 2 for any p, we can get the equation 
of Poo as a function of, p, to, n and (k)i 



pll-expi-ik^Poo)]* 



(32) 



The average number of cascading stages, (r), as a func- 
tion of p for different value of average degree is shown in 
Fig. 11. The numerical simulation results show that r 
increases sharply when p is near p c [20] . 
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FIG. 8: Loopless NON is composed of (a) ER networks, 
(b) RR networks and (c) SF networks. Plotted is P*, 
as a function of p for k = 5 (for ER and RR networks) 
and A = 2.3 for SF networks and several values of n. 
The results obtained using Eq. (27) for ER networks, 
Eq. (40) for RR networks and Eq. (51) for SF networks, 
agree well with simulations. 



IV. ANALYTICAL RESULTS FOR THE CASE OF 
NON OF n RR NETWORKS 

Next, we study the case for a tree- like NON of n RR 
networks. The degree of network i is fcj. The generating 
functions of network i are 



Gifi{x) = ^2 kiX 



and 



G ill (x)=x k *- 1 . 
Using Eqs. (2) and (3), we obtain, 

G ifi (f,x) = (fx + l-x) ki 



Gi,i(f,x) - (fx+l-x) 



ki-l 



(33) 
(34) 

(35) 
(36) 




FIG. 9: Loopless NON is composed of (a) ER networks, 
(b) RR networks and (c) SF networks. Plotted IS Poo 
a function of p for n = 5 for several values of k (ER and 
RR networks) and several values of m (SF networks for 
A = 2.3). The results obtained using Eq. (27) for ER 
networks, Eq. (45) for RR networks and Eq. (56) for SF 
networks, agree well with simulations. 




FIG. 10: For star-like network of 5 ER networks, the 
average convergence stage (r) as a function of p for dif- 
ferent k. In the simulation, N = 10 6 , and the simulation 
results are obtained by over 30 realizations. This feature 
enables to find accurate estimate for p c in simulations 
[20].) agree well with simulations. 



For a single network i we obtain gi{x) = 1 — Gi.o(fi, x), 
where fi satisfies the equation fi — Gi t i(fi,p). For loop- 
less NON of n networks, we can obtain 

n 

ti=p n [i -(./>;■ + (37) 
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FIG. 11: For loopless network of n ER networks, {k)\p c 
and (fc)i/u(oo) as function of the ratio {k)\/ {k)^ for n = 2 
(dashed), n — 3 (dotted), n = 4 (dashdot) and n = 5 
(solid) and for m = 1 (o), m = 2 (□), m = 3 (o) and 
m = 4 (>), where m denotes the number of individual 
networks whose average degree are the same (k)2 and 
average degree of the other n — m networks arc 
The results are obtained using Eqs.(29)-(31). 



where fj satisfies 
Thus we can obtain 

p 00 = l i 00 = P \{(i-f~). 

3 = 1 



(38) 



(39) 



We can obtain -Poo 

as a function of r by substituting r 

into Eq. (41), 

Poo ^ P {i-{^p^ i [(i-(^)*)^-i]+i}T- (45) 

The results to compare the simulation and theory are 
shown in Fig. 6(a). The numerical simulation results is 
shown in Fig. 8(b) and 9(b). 

For n ^> 1, Eq. (42) can be rewritten as 



y kn' 



(46) 



From Eq. (43) and Eq. (46) for the case when n » 1, 
we obtain [34] 



Pc = 



l-r 



(47) 



where /j satisfies 



where r satisfies Eq. (46). 

Since fc/(l — k) < for k > 1, it follows, in contrast 
to the ER case, in the RR NON case p c can never be 
greater or equal to 1. This shows that an ER NON is 
extremely more vulnerable compared to RR NON, due 
to the critical role played in the ER by singly connected 
nodes. 



V. ANALYTICAL RESULTS FOR THE CASE OF 
NON OF n SCALE-FREE NETWORKS 



h=u-i)p n (w/^+i]*'- 1 - (40) 

When all n networks have the same degree k i.e., fcj = k 
(i = 1,2, ...,n), we introduce a new variables r = f 1 *^ 
into Eq. (40), and the n equations are reduced to a single 
one 



Here we study the case of a tree-like NON composed 
of n scale-free(SF) networks. The generating function of 
each network is 



Glfi[x) ~ (M + 1)1-* -m*-* ' (48) 



r = ( r *-i - i) p (i - r fc )«-! + l 



(41) 



which can be solved graphically for any p. The critical 
case corresponds to the tangential condition. Thus, we 
obtain that the value of r satisfies a transcendental equa- 
tion 



1 = (l-r )r k 2 [-- >— + - 

1 — r K 1 — r 



k-l J 



(42) 



Solving r from Eq. (42), we can obtain the critical value 
of p c and the the value of at p c as 



and 



cM- ^KHil^-f-V- 1 (49) 



G il0 (/,a;) = G i ,o(l + a;/-a;) ) (50) 



Ga(/,^Gu(l + x/-4 (51) 



r - 1 



(j.k—1 _ 1)(1 — r k\n— 1 ' 

The numerical solutions are shown in Fig. 7(b). 

Poo =p c (l-r fc )". 



(43) 



(44) 



l-G ifi {l + xf -x). 



fi = G l . 1 (l + xf -x). 



(52) 



(53) 



Substituting Eq. (50) - (53) into Eqs. (8)-(10), wc 
obtain 
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T-r r Eg [(fc + 1) 1 -" - fc^Ki + fjxj xtf 1 

x Pill 1 (M + l^-m 1 -^ J ' 

(54) 



/i 
and 



E^Kfc + i) 1 -*-* 1 -*]* 



(55) 



71 

n{i- 



l-A 



„1-A 



](! + /: 



(M + l)!- A 



l-A 



(56) 

Figs. 8(c) and 9(c) show the solutions for for sev- 
eral values of n and m respectively. 



VI. THE CASE WHEN n ->• oo 

The behavior of the networks for n — > oo depends only 
on P(k), for small k. 

If P(l) + P(0) > 0, then no matter how large is (k) 
and what is the degree distribution for the rest of k, the 
networks completely collapse for large enough n (Poo = 

o). 

If P(l) + P(0) = the networks survive (P^ > 0) for 
any n and large enough p < 1. Indeed, if P(2) > but 
P(l) + P(0) = 0, 



i( 



1 



J_) 

n v 2P(2) 



(57) 



which corresponds to the maximum of the r.h.s of Eq. 
(12) and p c — > 1, P^ — >• 1 can be found from Eqs. (12) 
and (13). Thus, p c can not be greater than 1, meaning 
that for all n values, the NON is stable. 



Assuming £o 1 P(k) = and P(£) > then, 

z c ~ [l/iniPii))} 1 ^- 1 ^. (58) 

and the way p c — > 1 and Px, — > 1 can be found from Eqs. 
(12) and (13) to be. 

= 1-Ci/n 1 /^- 1 ) andp c = 1 - Cs/n 1 /^" 1 ), where 
C\ > and C2 > are constants, that can be easily 
found from Eqs. (12) and (13). 

For P(2) > 0, C 2 - 2P(2)/(fc) 2 and Ci = (1 - 
2P(2)/(fc)) 2 /(2P(2)) and for 

P(2) = 0, C 2 = (£ - l)/(£P(£)Y^ e -^ and d = 
1/(£P(£))] 1 M- 1 '>. Thus, when P(0) + P(1) = 0, the NON 
is stable for all n (p c (n) < 1) and a condition for a mini- 
mal k(n), such as in Eq. (26) does not exist. 

VII. CONCLUSION 

In summary, we have developed a framework, Eqs. (9)- 
(10), for studying percolation of NON from which we de- 
rived an exact analytical law, Eqs. (27) [for ER networks] 
and (45) [for RR networks] , for percolation in the case of 
a network of n coupled networks. In particular for any 
n > 2, cascades of failures naturally appear and the phase 
transition becomes first order transition compared to a 
second order transition in the classical percolation of a 
single network. These findings show that the percolation 
theory of a single network is a limiting case of a more 
general case of percolation of interdependent networks. 
Due to cascading failures which increase with n, vulnera- 
bility significantly increases with n. We also find that for 
any tree-like network of networks the critical percolation 
threshold and the mutual giant component depend only 
on the number of networks and not on the topology (see 
Fig. 2). We discuss the case for n coupled ER networks, 
RR networks and SF networks. We find that there ex- 
ist the minimum k to make the NON survives, but no 
parameter for the RR and SF networks. 
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introduce a new variable 



7 ■ 



nfc 



(61) 



so, f c — e~ 7 . When n — > oo, f c — > and 7 — > 00. So 
studying the case n — > 00 is the same as studying the 
case 7 — > 00. 

Substitute Eq. (20) to Eq. (19), we obtain 



1-e" 7 
n = . 

7e-f 

We study n as a function of 7 when 7 — > 00, 



lim n = — 

7^00 7 



Inn = 7 — In 7. 



(62) 



(63) 



(64) 



Substituting Eq. (21) to Eq. (16) and consider the case 
when 7 — > 00, we obtain 



fcmin — 



(1 - e-T)™ : 



lim k n 

7—5-00 



(65) 



(66) 



VIII. APPENDIX 



A. The case when n -¥ oc for ER NON 



Substituting Eq. (23a) to Eq. (22b) and consider the 
case when 7 — > 00, we obtain 



Eq. (13) can be written as 

fc-l 



n 



fc In f c 



Then we can get, 



dn In f c + 1 - f c 



dfc (/cln/c) 



2 ' 



(59) 



(60) 



When n = l,f c = 1, and 4£ = 0; when n > l,/ c < 1, 



and ^- < 0. 



So f c is a decreasing function of n when n > 1. We 



fcmin = e Tlnn + e 7 In fc min + -e -< , (67) 

7 



lim fc m in = ln(nlnn) + £, 

7— >oo 



(68) 



where ^ = 0(ln(lnn)). We can also obtain that when 
n — > 00, p c satisfies 



Pc = -ln(nlnn) +C, 



(69) 



where £ = O (In (Inn)). This result is corresponds to large 
n values in Fig. 7(a). 
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